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A MODEL FOR UNDERSTANDING NUMERICAL STABILITY 


POLKMAR BORNEMANN* 

Abstract. We present a model of roundoff error analysis that combines simplicity with predic¬ 
tive power. Though not considering all sources of roundoff within an algorithm, the model is related 
to a recursive roundoff error analysis and therefore capable of correctly predicting stability or in¬ 
stability of an algorithm. By means of nontrivial examples, such as the componentwise backward 
stability analysis of Gaussian elimination with a single iterative refinement step, we demonstrate 
that the model even yields quantitative backward error bounds that show all the known problem- 
dependent terms (with the exception of dimension-dependent constants, which are the weak spot of 
any a priori analysis). The model can serve as a convenient tool for teaching or as a heuristic device 
to discover stability results before entering a further, detailed analysis. 
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1. Introduction. An algorithm for the numerical evaluation of a complicated 
function f is just a decomposition into simple intermediate steps, such as arithmetic 
operations, elementary transcendental functions, or well-behaved and well-understood 
library algorithms (e.g., matrix multiplication): 

f = gi o 92 o • • • o gk. 

In the realm of floating point arithmetic, each of these intermediate steps is contam¬ 
inated by roundoff and hence contributes to the final perturbation of the result in a 
twofold fashion: first, by generating roundoff error itself, and second, by propagating 
the roundoff errors of previous steps. Since the early days of numerical computing 
there has been much progress in clarifying the underlying structure and organiz¬ 
ing the results in a concise, easily interpreted form. However, a detailed analysis 
(Higham 2002) is still often quite involved and remains a battle-field for experts, too 
tedious to teach and explain in detail beyond the most trivial cases in a beginner’s 
course on numerical analysis. The instructor typically chooses between two options: 
skipping the nontrivial results (such as stability of Gaussian elimination) entirely, or 
just stating the results without proof. Either choice is unsatisfactory for good stu¬ 
dents since they cannot develop an understanding of the mathematical structure and 
reasons. 

We will demonstrate in this paper, that the overall behavior of an algorithm can very 
often be well understood by analyzing a simplified model of the sources of roundoff 
error. As in the natural sciences such a model has to balance simplicity with predictive 
power. If such a simple model basically leads to the same predictions, qualitatively 
and perhaps even quantitatively, as a full-fledged a priori roundoff error analysis, 
we may rightly claim to have contributed to the understanding of the algorithm’s 
behavior. In fact, all the estimates of our model analysis that we present in this paper 
will give the same estimates as a detailed a priori analysis—with the only exception of 
the dimension-dependent constants, which are, however, anyway the weak spot, and 
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therefore least important part, of any roundoff error analysis (Higham 2002, p. 65). 
In particular, with just a few lines of simple calculations we will obtain the nontrivial 
results on the norm- and componentwise backward stability of Gaussian elimination 
ranging from the early work (Wilkinson 1963) to the analysis of iterative refinement 
(Skeel 1980). 

In addition to being a convenient (and to the experience of the author also successful) 
tool in teaching, our model might serve as a heuristic device in discovering the 
structure of a stability result—before one enters, in a second step, taking advantage 
of the obtained knowledge, a fully detailed roundoff error analysis. 

The Model. The roundoff error analysis that we propose is based on the obser¬ 
vation that in many if not most cases a critical intermediate step can be identified 
that leads to a natural decomposition 

f = gi o • • • o gj o gj + i o • • • o gk = ho g 
=h =g 

into just two fundamental steps. Now, the model is based on the simplifying as¬ 
sumption that roundoff error just affects the single intermediate result—after being 
output by g, before being input to h. That is, we analyze the error of the realization 
map 

f = h o ft o g. 

Here ft: denotes componentwise rounding, subject to the standard model 

of floating point arithmetic 

|fl(x) — x| < u • |x|, 

where u denotes the unit roundoff of the arithmetic (u « 1.11 X 10-1*5 for IEEE double 
precision) and F the fioating point numbers. We understand jx| componentwise for 
vectors and matrices. 

Outline of the Paper. In |2lwe analyze the backward stability of the realization 
map f, which turns out to be determined by the condition number of g-'. We will 
specify the relation of the model to a complete analysis. In fact, if the model is 
unstable the same has to be expected for the real situation. On the other hand, if g 
and h are realized by backward stable algorithms, the resulting algorithm for f would 
inherit the stability of the model. This helps to understand the success of our model 
and suggests a recursive approach to a full roundoff error analysis. 

The rest of the paper studies some algorithms for the solution of a linear system 
Ax = b. In we recall some classic expressions for the backward error of linear 
systems that are the point of departure for the simple estimates to follow. In we 
study the naive algorithm, that is, multiplication with A-^, and show its instability 
for badly conditioned matrices. In 33 we study the normwise backward error of 
Gaussian elimination and obtain the classic result (Wilkinson 1963). In IH^e get the 
result (Skeel 1979) on the componentwise backward error of Gaussian elimination, 
correctly predicting the influence of the scaling of the system. Finally, in 0 we show 
how to discover within the frame of our model the result (Skeel 1980) that a single 
step of iterative refinement implies componentwise backward stability of Gaussian 
elimination. 
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2. Backward Stability. The main result of a qualitative study of our model 
can be summarized as follows: 

Backwards stability requires that is well conditioned. 

In fact, backward stability analysis requires the result of the algorithm for input x, 
that is f(x) in our model, to be represented as the exact solution to perturbed data: 
f(x) = f(x + Ax). Writing w = g(x) for short, we have 

f(x) = h(fl(w)) = h(w + Aw), |Aw| = |fl(w) — w| < u • |w|. 

Assuming g to be invertible, we propagate Aw backwards to obtain an estimate for 
Ax: 


x + Ax = g ^(w + Aw), |Ax| < Kg-i u • |x| + O(u^), 

where Kg-i denotes the (componentwise) relative condition number of g^^ at w. 
Hence, the backbard error is bounded by the unit roundoff amplified by Kg-i. 

2.1. Examples. 

A. Consider the evaluation of f(x) = log^(1 + x) for x « 0. A direct implemen¬ 
tation of the defining formula corresponds to the decomposition 

f : X 1-^ w = 1 -|- X log^ (w). 

Now, because w = 1 -|-x « 1 the inverse function g^^ : w i—> x = w— 1 is a subtraction 
in the cancelation regime, thus badly conditioned. Hence, we predict instability of the 
formula, which simple examples confirm. In fact, the bad conditioning of g^^ reflects 
the loss of information in g: we have fl(g(x)) = fl(1 -1-x) = 1 as soon as x is smaller 
than the resolution of the machine arithmetic. In general, well-conditioning of g^^, 
however, requires that the input x is accurately reconstructable from the intermediate 
result w = g(x). 

B. The solution x G of a linear system of equations Ax = b with a nonsingular 

A G can formally be respresented as x = A^^ • b. This suggest the naive 

algorithm corresponding to the decomposition 

f : A A-^ X = A-’ • b. 

Now, g^^ : A^^ h—> A is just g again, its condition is (in the normwise case) the 
condition number of the matrix A. Thus, we expect the algorithm to be unstable for 
certain badly conditioned matrices. Examples that display such instability will be 
given in where we extend our analysis to a more quantitative setting. 

C. On the other hand, the solution of the linear system Ax = b by Gaussian 
elimination corresponds to the decomposition 

f : A (L, U) X. 

Here, g represents the LU-factorization step, whereas h represents the substitution 
steps. Now, the inverse of g, that is 
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is just matrix multiplication. Its condition number can be estimated by 

|L|-|U||| 


Kg-l < 2- 


iA|| 


which is, as will be discussed in more detail in iH sufficient to explain the instabilities 
to be observed for Gaussian elimination with or without partial pivoting. 

2.2. Relation of the Model to a Complete Analysis. In fact, the condition 
number of turns out to be relevant for a full roundoff error analysis, too. Here, 
we would recursively define the realization of f = K o g by 

f = ho g, 

starting with the backward stable realization of the arithmetic operations and basic 
elementary functions. (Of course, in general we cannot assume that in each step of 
this recursion the g-part of the decomposition is invertible. However, it is possible to 
give a reasonable definition of Kg-i even if g is many-to-one.) 

With |AxJ| denoting maximum componentwise relative error^ we define the smallest 
number |3f > 0 such that 

f(x] = f(x + Ax), |Ax| < |3f • u+O(u^) 

as the stability indicator of f. Backwards stability requires (3f to be not too large. 
Lemma 2.1. For g invertible there holds the recursive estimate 


Pf ^ Pg + Kg- 

Proof. The stability indicator of h gives 

f(x) =h(g(x)) =h(g(x) +Aw), 


Ph- 


|Aw] < Ph -uT O(u^). 


( 2 . 1 ) 


The stability indicator of g and the relative condition number of g ^ allow for the 
estimates 


g(x) + Aw = g(x + Axi) + Aw, 
= g(x + Axi + Ax 2 ), 


|Axi] < Pg -uT O(u^), 

[Axil < Kg-, • |Aw] + 0(|Awf 


Since Axi and Ax 2 are both perturbations of the same quantity x there holds the 
triangle inequality for relative errors. 


Ax = Axi + Ax 2 , 


|Ax]| < |Axi] + |Ax2] < (Pg + Kg-i • Ph)u+ O(u^), 


and we get the assertion. □ 

Thus, we may complement the maxim from the beginning of this section by the 
following rule: 

If g^^ is well-conditioned, backward stable realizations of g and 
h induce a backward stable realization o/f. 


^That is, for the perturbation Ax S 1 
with the convention that 0/0 = 0. 


of a quantity x £ we have [Ax] = maxj = i Axj |/|Xj I 
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Summarizing, the logical status of the proposed model is a follows. If the model 
predicts instability, we can expect instability in reality—independently of how g and 
h are realized in practice. Most probably, in examples that realize the worst case 
scenario of the condition number bound, there will be instability even if g and h 
were calculated exactly; a fact, which certainly shakes our faith in the algorithm. On 
the other hand, if the model predicts stability, the actual stability of the algorithm 
depends on how g and h are realized algorithmically. In the framework of back¬ 
ward stability, stability of the realization of g and h implies stability of the resulting 
algorithm for f. 

Example. Let us illustrate these points by reconsidering the example of 92.II A. 
Here, we decompose f(x) = log'll! -|- x), x « 0, differently into 

f : X 1-^ w = log(l -|- x) w^. 


Now, the critical map g^^ : w i—> e™ — 1 has relative condition number Kg-i ft! 1 for 
w ft 0. The model alone would therefore predict numerical stability. On the other 
hand, the full, recursive analysis has to take the actual algorithms for g and h into 
account. Step h, as a multiplication in IEEE arithmetic, is certainly backward stable. 
However, the status of g is far less clear. If its realization is chosen to be based on 
the decomposition gtxi—>z=l-|-xi—> log(z), then an analysis similar to 32.H A 
reveals instability. Otherwise, if g is realized, for instance, by using Kahan’s stable 
algorithm as implemented in Matlab’s loglp command, the resulting algorithm for f 
is stable, too. 

Hence, the choice of the decomposition will critically determine the success or failure 
of the model. In general, making a conclusive choice will depend on the user’s ex¬ 
perience or luck. However, we will show in the rest of the paper, that quite natural 
such decompositions occur in the analysis of the stability of Gaussian elimination. 

3. The Backward Error of Linear Systems. To prepare for a more quanti¬ 
tative analysis of algorithms for the solution of linear systems of equations Ax = b 
we recall the concept of the backward error of an output vector x S Normwise 
analysis considers^ 

ri = min | : (A E)x = b 1 , 

egR'^x™ |jA|| J 

whereas componentwise analysis studies 

cu = min < max 7 — 7 ^ : (A -|- E)x = b 

EgRmxm 1 ^ ij |A|7j 


The classic results (Rigal and Caches 1967) and (Oettli and Prager 1964) show that 
q and cu can be calculated from the data of the linear system and the output vector 
X by means of the following simple formulas: 


h = 


l|r|| 

l|A|!-|!x||’ 


CU = max 

j = ^ :m 


Ip 

(|A!-|x|)j- 


(3.1) 


^Throughout the paper we deal with monotone vector norms like the 1-, 2-, or oo-norm, and the 
induced matrix norms. 
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Here, r = b — Ax denotes the residual of x. These formulae, which have very short 
and straightforward proofs (Higham 2002, pp. 120/122), are also valuable for the a 
posteriori assessment of computed solutions. We will use them as a convenient point 
of departure for a quantitative analysis in the frame of our proposed model. 

4. Model Analysis of the Naive Algorithm for Linear Systems. As dis¬ 
cussed in il2.1l the naive algorithm for the solution of a linear system is given by the 
decomposition 

f : A B = X = B • b. 


Our model analyzes how roundoff in B affects the solution x and its backward error: 


f: Ai 


B = A- 


B = B 


AB 


-!^x = B-b. 


The perturbation |AB| < u • |B| induces, by propagating backwards through , an 
equivalent perturbation A = A- 1 - AA = (B) of the input matrix. By construction, 

we have Ax = b, 

(A-F AA)(A^^ AB) = I, i.e., AA =-A • AB • A - AA • AB • A, 
and therefore the componentwise estimate 

|AA • x| < |A| • |A-^ I • |Ax| u -f 0(u2). 


Since t = b — Ax = AA • x and x = x + 0(u), we get by 

-’l-IAxll' 


||AA.x|| ^ IIIAI-IA- 

ri = < 


■u-h 0(u ) y(A,x)u-|- 0(u ). 


l|A||-||x|| - ||A||.||x|| 

To relate with better known quantities, we may further estimate 

y(A,x) < li |A| • lA^'l II = cond(A^’'], 


(4.1) 


in agreement with our qualitative analysis of 82.1I B. Thus, instability in the normwise 
concept appears to be only possible for badly conditioned matrices. 

4.1. Examples.^ 

A. A notoriously badly conditioned matrix is the famous Hilbert matrix Hra for 
larger dimensions m. In Matlab there is the command invhilb that supplies 
and allows to implement the naive algorithm:'* 


>> m = 20; A = hilb(m); B = invhilb(m); b = ones(m,l); x = B+b; 
>> eta = norm(b - A*x,inf)/norm(A,inf)/norm(x,inf) 


eta = 1.2787e-005 


®If not explicitly stated otherwise, all the examples in this paper use the norm || ■ ||oo. 

^Here, and in the examples to follow, we have cross-checked the actually calculated backward 
errors with higher precision arithmetic. The first digits were always correct, so that the conclusions 
we draw are not affected by roundoff errors in the computed residuals. 
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Thus, the naive algorithm is unstable as predicted by the a priori bound which 

turns out to be 


r| = 1.27--- X 10^^ <y(A,x) x 10^"^; 

a fairly good prediction indeed. On the other hand, we have to be careful to base a 
prediction on coarser upper bounds that were introduced for the ease of interpretation: 
the condition number yields 

T] < cond( ] • u = 6.63 • • • x 10^ \ 

which gives too pessimistic a picture of the actual backward error. 

B. The following example (Skeel 1979, p. 509) shows that the naive algorithm 
can be stable for some badly conditioned matrices: 


/I 1 -1 -1\ 




(w 

10 0 1 

■U _ 

2 


e ^ 

0 0 e 0 

) L) — 

1 

) ^ — 


Vo e 0 0 j 




^ 1 / 


This matrix fulfills 

cond(A) = 4, cond(A^^) = 2 + 4e^'. 

However, numerical experiments with various small 0 < e <C 1 exhibit very small 
backward errors of about the size of the unit roundoff. This is fully refiected by our 
model analysis, since 

y(A,x) = 1 + I « 1. 

5. Model Analysis of Gaussian Elimination: The Normwise Case. As 
discussed in 32.1I C the solution of a linear system Ax = b by Gaussian Elimination 
corresponds to the decomposition 

f : Ah-?^ (L,U) h-!^x. 

In the model roundoff affects only the intermediate result, the LU-factorization, by 

f : A (L, U) (L, U) = (L + AL, U + AU) x. 

Here, the perturbations |ALj < u • |L|, |AU| < u • |U| induce, by propagating through 
the inverse of g (that is, matrix multiplication), an equivalent perturbation of the 
input matrix 

A + AA = LU= (L + AL) • (U + AU), i.e., AA = AL • U + L • AU + AL • AU. 
This way we obtain the componentwise estimate 


|AA| <2|L||U|.u,, 


u* = u + u^/2. 


(5.1) 
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Because of r = b — Ax = AA • x we get by 


|||AA|.|x||| ^ J||L|.|U|.|x||| ^ J||L|.|U||| ^ 

|lA||-||x|l ||A||-||x|i ||A|| . 27 ( 1 ,U)u,, (5.2) 

in agreement with our qualitative analysis of il2.1I C. If we restrict ourselves to mono¬ 
tone matrix norms, we can further estimate the growth factor y(L,U) by using 
U = L-i A 




|A|| 


All 


= cond(L ^I 


Thus, an instability of Gaussian elimination in the normwise case requires a badly 
conditioned L-factor of the matrix A. 

5.1. Examples. 

A. It is well known that Gaussian elimination without pivoting is bound to be 
unstable for small pivot elements. An example is given by 


A = 


e 1 
1 1 


L = 


1 0 
1 


U = 


e 1 

0 1-e-i 


For e = u, b = (1,0)^ a numerical experiment yields® x = (—2,1); the exact solution, 
however, would be x = (—1,1 )^. The backward error turns out to be p = ;j. On the 
other hand we have 


Y(L,U)=e \ cond(L ’)=1-|-2e \ 

which, by H5.2|l . gives the fairly good prediction p < e ^ • u = 1. 

B. Gaussian elimination with partial pivoting yields an L-factor that satisfies 
|L| < 1 componentwise. This can be used (Higham 2002, p. 143) to show that 

y(L,U) < cond(L-M < 2’^-!, 


which proves that the growth factor remains bounded for fixed dimension m. How¬ 
ever, the upper bound on cond(L^^) is attained for Wilkinson’s famous matrix 

\ 

1 1 / 

Numerical experiments quickly exhibit very large backward errors: 


/ 1 


A = 




L = 


V -1 ••• -1 1 / 


( 1 
-1 

V -1 


>> m = 53; A = eye(m)-tril(ones(in) ,-l) ; A(:,m) = 1; 
>> rand(’seed’,42); b = rand(m,l); x = A\b; 

>> eta = norm(b-A*x,inf)/norm(A,inf)/norin(x,inf) 


eta = 3.2342e-003 


®We write a = bifa — b^u. 
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Our analysis yields a fairly good prediction, 

11 =3.23--- X 10-3 < 27 (L,U) -u* =3.77--- x IQ-^. 


C. For symmetric positive definite matrices, the solution of the linear system 
Ax = b by Cholesky factorization corresponds to the decomposition 

f: Ah-^Lhii^x 

with A = L • A perturbation L = L + AL of the intermediate result by roundoff. 


|AL| <u-|L|, 

induces, as for l|5.2|l . the backward error (with resprect to the norm 11-112) 

"|L|-|LT|||2 


il<2- 


llAli 


■u, =2 y(L, 


Since || |L| II 2 < •\/Ta||L|j 2 for any m x m matrix, we infer (Higham 2002, p. 198) 


Y(L,LT)< 


|l-lll2|||L’'l 

IIAIU 


< m- 


II2||L^||2 

Illtil 


= m. 


Hence, we have 


U < 2mu*. 


which hints to the perfect normwise backward stability of the Cholesky method. 

6. Model Analysis of Gaussian Elimination: The Componentwise Case. 
The matrix estimate JOl immediately yields an estimate of the componentwise 
backward error, 


tu 


max 


|AA ■ x|j 
(|A| • |x|)j 


< max 


(|AA| ■ |x|)j 
(|A|-|x|)j 


< 2 max 


(|L| • |U| • |x|)j 
(|A|-|x|)j 




^ maxj(|L| ■ |U| ■ |x|)j ^ || |L| ■ |U| ■ |x| ||^ maxjdAj ■ |x|)j 

minj(|A| • |x|)j) II !A| • |x| ||^ miuj (|A| • |x|)j 


which by U = A, that is |U| < |L^'| • |A|, induces (Skeel 1979, Thm. 4.4) 

cu < 2 cond(L-^) cr(A, x) u*. (6-2) 

As our derivation shows, this is not necessarily the best possible concise bound, but 
it allows for the easy comparison with the normwise bound (with respect to ll'Hoo) 

r) < 2cond(L-^) u*. 


We see that the componentwise bound just differs by the additional factor cr(A, x) > 1. 
This factor measures the quality of the scaling of the linear system with respect to x 
and predicts an instability for badly scaled systems. 
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6.1. Examples. 

A. We return to the example of 94.1I B. The growth factor and the scaling are 
given by 


cond(L ^)=3+4e, a(A,x')—2 + 2e \ 

Experimentally, for e = Gaussian elimination yields (partial pivoting is not 

used here because of |L| < 1) 

ri =2.84--- X 10-^^ <2cond{L-^)u* =6.66--- x 10^^^ 

On the other hand, the componentwise backward error satisfies 

cu = 0.499 • • • < 2cond(L^’ (o’!A, x)u* = 13.3 • • • . 


Thus, the model analysis helps to understand the actual behavior of the two error 
concepts. In particular, we see that scaling can be an issue for Gaussian elimination 
with partial pivoting if analyzed componentwise. 

B. There are matrices, for which the upper bound l|6.2|l turns out to be too 
coarse. As an example, we consider totally positive matrices A such as the Hilbert 
matrix of 94. II A or matrices that appear in spline interpolation. These matrices factor 
with L > 0 and U > 0. Thus, we best stay with the following intermediate step in 
the chain of estimates 


cu < 2 max 

j 


(|L| ■ |U| • |x|)j 
(|A|-|x|)j 


u*. 


Here, we obviously have |L| • |U| = |A| and we can therefore directly infer the perfect 
stability estimate (de Boor and Pinkus 1977) 


cu < 2u*. 


7. Model Analysis of a Single Iterative Refinement Step. In this final 
section we will apply the model analysis to the understanding of the results (Skeel 
1980) on iterative refinement of Gaussian elimination. We recall that the iterative 
refinement of a calculated solution x of a linear system Ax = b consists of three 
steps: compute the residual to = b — Ax, solve Aw = Tq for a calculated correction 
w (reusing the LU-decomposition of A), update y — x + iv. If there were no roundoff 
errors in the refinement steps (that is, w = w), we would obtain y = x, the exact 
solution. 

In the previous two sections, the model analysis of Gaussian elimination allowed for 
roundoff errors just in the L- and U-factors of A yielding some equivalent perturbation 
of that matrix. Because of the reuse of these factors in the iterative refinement step, 
we reasonably assume that both Gaussian elimination steps, that is, those leading to x 
and w, are affected by roundoff through a single perturbation A = A +AA satisfying 
the estimate in. This way, the result y of the iterative refinement is given by 


y — x + w 


(A + AA)w = To = b — Ax = AA • x. 
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The residual after this step is ti = b — Ay = Tq — Aw = AAw, and therefore 
Aw = To — AAw = AA(x — w) = AA(y — 2w), w = A^^ AAy — 2A^' AAw. 
Hence we have 


|AAw| < |AAj|A-^||AA||y| + 2 |AA| lA^^ | |AAw|, 
which by (jOJ, that is |AA| < 2 |Lj |U| u, < 2 |L| | |A| u*, implies 

||AAw||oo < 4cond^(L^’') cond(A^’') • 11 |Ai |y| 11^ uj 

+ 4cond(L^’) cond(A^^ )u* ||AAw||^ . 

If 4 cond(L^') cond(A^^) u* < 1 we can solve for ||AA w|joo and get—as in the deriva¬ 
tion of l|6.1|l —the following upper bound of the backward error of y: 


CU] 


|ri Ij |AA • w|j maxj|AA-w|j 

T (|A| • |y|)j ^ “f"" (|A| • |y|)j " minj (|A| • |y|)j 

^ 4cond^(L~^) cond(A~^) g(A,y)uH. 

l||A||y|||;^ “ 1-4cond(L-’) cond(A-’)u* 


Because of cond(L“ -’) > 1, o-lA.y) > 1, the premise is in particular satisfied if 


8cond^(L ^) cond(A ’) crfAj-y)u* < 1, (7.2) 


for which we obtain from jni the simple perfect bound cui < u*. Except for a 
constant depending on the dimension m this is exactly the result (Higham 2002, 
p. 239) of an elaborate analysis that takes all the details of roundoff error rigorously 
into account.® 

Summarizing our analysis predicts: As long as the linear system is not too badly 
conditioned (cond(A^^) is not too large) and not too badly scaled (cr(A,'y) is not too 
large), and Gaussian elimination is not too unstable (cond(L^') is not too large), one 
step of iterative refinement implies componentwise backward stability. 

7.1. Example. We consider the example (Skeel 1979, p. 500) 


A = 


2 

2e 

2e 




x = 



of a well conditioned (for this particular right hand side b), but badly scaled linear 
system. Because of 

(2 "I Q 

cond(A^^) =-|-0(1), o'(A,x) =+ 0(1), cond(L^^) = --|-0(e), 

®The catch, of course, is that without doing the full analysis we would not know if we had 
really determined the full bound. However, the point of this paper is a better understanding of the 
underlying mathematical structure. If, by neglecting many details, we come to predict the same 
bounds with much less effort we seem to have put the focus on the right spot. 
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Fig. 7.1. Backward errors ujq and cu] vs. e 


condition reads as 

')C/1 

1 > 8 cond^(L^^) cond(A^’) crtA,^) u* = —+ 0(e^^), 

5 

that is, one step of iterative refinement is predicted to imply stability as long as e 
remains larger than about the square root of the unit roundoff, 

e > ^V^+0{u^) « 7.5 X 10-*. 

5 

In fact, the upper bound l(0|l predicts that the componentwise backward error cuo of x 
behaves like cuq = O(e^^u); whereas the upper bound jni predicts cui = 0(e 
for the first refinement step y. All this can perfectly be observed in an actual numerical 
experiment, see Figure [TTTl 
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